/*=========================================================================
 *
 *  Copyright Insight Software Consortium
 *
 *  Licensed under the Apache License, Version 2.0 (the "License");
 *  you may not use this file except in compliance with the License.
 *  You may obtain a copy of the License at
 *
 *         http://www.apache.org/licenses/LICENSE-2.0.txt
 *
 *  Unless required by applicable law or agreed to in writing, software
 *  distributed under the License is distributed on an "AS IS" BASIS,
 *  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
 *  See the License for the specific language governing permissions and
 *  limitations under the License.
 *
 *=========================================================================*/

// Software Guide : BeginLatex
//
// This example illustrates the use of the \doxygen{DeformableMesh3DFilter}.
// An initial mesh is created using the \doxygen{SphereMeshSource} filter.
//
// \index{Deformable Models}
// \index{DeformableMesh3DFilter}
// \index{SphereMeshSource}
//
// Software Guide : EndLatex

// Software Guide : BeginCodeSnippet

#include "itkDeformableMesh3DFilter.h"

#include "itkGradientRecursiveGaussianImageFilter.h"
#include "itkGradientMagnitudeRecursiveGaussianImageFilter.h"


#include "itkSphereMeshSource.h"

#include "itkImageFileReader.h"

int main(int argc, char * argv [] )
{

  if( argc < 7 )
    {
    std::cerr << "Usage: DeformableModel2 inputImage sigma " << std::endl;
    std::cerr << " numberOfIterations timeStep externalFoceScale stiffness" << std::endl;
    return 1;
    }
  const char *inputFileName      =       argv[1];
  const float sigma              = atof( argv[2] );  // Suggested value = 5 * pixel spacing
  const int   numberOfIterations = atoi( argv[3] );  // Suggested value = 100
  const float timeStep           = atof( argv[4] );  // Suggested value = 0.1
  const float externalForceScale = atof( argv[5] );   // Suggested value = 10  (linked to stiffness)
  const float stiffness          = atof( argv[6] );   // Suggested value = 0.1 (linked to force scale)


  typedef double                            MeshPixelType;
  typedef itk::Mesh< MeshPixelType >        MeshType;


  unsigned const int Dimension = 3;

  typedef   float                               PixelType;
  typedef itk::Image<PixelType, Dimension>      ImageType;


  typedef itk::GradientRecursiveGaussianImageFilter<
                                        ImageType
                                           > GradientFilterType;

  typedef itk::GradientMagnitudeRecursiveGaussianImageFilter<
                                        ImageType,ImageType>
                                                      GradientMagnitudeFilterType;


  typedef itk::DeformableMesh3DFilter<MeshType,MeshType>  DeformableFilterType;


  typedef itk::ImageFileReader< ImageType       >  ReaderType;

  ReaderType::Pointer       imageReader   =  ReaderType::New();

  imageReader->SetFileName( inputFileName );

  GradientMagnitudeFilterType::Pointer  gradientMagnitudeFilter
                      = GradientMagnitudeFilterType::New();

  try
    {
    imageReader->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Problem found while reading file " << inputFileName << std::endl;
    std::cerr << "Exception Caught !" << std::endl;
    std::cerr << excep << std::endl;
    }


  ImageType::ConstPointer inputImage = imageReader->GetOutput();

  gradientMagnitudeFilter->SetInput( inputImage );


  gradientMagnitudeFilter->SetSigma( sigma );


  GradientFilterType::Pointer gradientMapFilter = GradientFilterType::New();


  gradientMapFilter->SetInput( gradientMagnitudeFilter->GetOutput());
  gradientMapFilter->SetSigma( sigma );

  gradientMapFilter->Update();



  DeformableFilterType::Pointer deformableModelFilter =
                                     DeformableFilterType::New();



  typedef itk::SphereMeshSource< MeshType >        MeshSourceType;

  MeshSourceType::Pointer meshSource = MeshSourceType::New();



  // Set the initial sphere in the center of the image
  //
  const ImageType::SpacingType spacing = inputImage->GetSpacing();
  ImageType::PointType         origin  = inputImage->GetOrigin();
  ImageType::SizeType          size    = inputImage->GetBufferedRegion().GetSize();


  MeshType::PointType center;
  center[0] = origin[0] + spacing[0] * size[0] / 2.0;
  center[1] = origin[1] + spacing[1] * size[1] / 2.0;
  center[2] = origin[2] + spacing[2] * size[2] / 2.0;
  meshSource->SetCenter( center );

  MeshType::PointType radius;
  radius[0] = spacing[0] * size[0] / 4.0;
  radius[1] = spacing[1] * size[1] / 4.0;
  radius[2] = spacing[2] * size[2] / 4.0;
  meshSource->SetScale( radius );


  meshSource->SetResolutionX( 50 );
  meshSource->SetResolutionY( 50 );
  meshSource->Update();


  deformableModelFilter->SetInput(    meshSource->GetOutput()        );

  deformableModelFilter->SetGradient( gradientMapFilter->GetOutput() );



  typedef itk::CovariantVector<double, 2>           StiffnessType;

  StiffnessType stiffnessVector;
  stiffnessVector[0] = stiffness;
  stiffnessVector[1] = stiffness;



  deformableModelFilter->SetTimeStep( timeStep );
  deformableModelFilter->SetStiffness( stiffnessVector );
  deformableModelFilter->SetStepThreshold( numberOfIterations );
  deformableModelFilter->SetGradientMagnitude( externalForceScale );


  try
    {
    deformableModelFilter->Update();
    }
  catch( itk::ExceptionObject & excep )
    {
    std::cerr << "Exception Caught !" << std::endl;
    std::cerr << excep << std::endl;
    }

  return 0;
}


// Software Guide : EndCodeSnippet
